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Abstract 

The  Cahn-Hilliard  equation  involves  fourth-order  spatial  derivatives.  Finite  ele¬ 
ment  solutions  are  not  common  because  primal  variational  formulations  of  fourth- 
order  operators  are  well  defined  and  integrable  only  if  the  finite  element  basis  func¬ 
tions  are  piecewise  smooth  and  globally  C  ^continuous.  There  are  a  very  limited 
number  of  two-dimensional  finite  elements  possessing  (^-continuity  applicable  to 
complex  geometries,  but  none  in  three-dimensions.  We  propose  Isogeometric  Anal¬ 
ysis  as  a  technology  that  possesses  a  unique  combination  of  attributes  for  complex 
problems  involving  higher-order  differential  operators,  namely,  higher-order  accu¬ 
racy,  robustness,  two-  and  three-dimensional  geometric  flexibility,  compact  support, 
and,  most  importantly,  the  possibility  of  C1  and  higher-order  continuity.  A  NURBS- 
based  variational  formulation  for  the  Cahn-Hilliard  equation  was  tested  on  two-  and 
three-dimensional  problems.  We  present  steady  state  solutions  in  two-dimensions 
and,  for  the  first  time,  in  three-dimensions.  To  achieve  these  results  an  adaptive 
time-stepping  method  is  introduced.  We  also  present  a  technique  for  desensitiz¬ 
ing  calculations  to  dependence  on  mesh  refinement.  This  enables  the  calculation 
of  topologically  correct  solutions  on  coarse  meshes,  opening  the  way  to  practical 
engineering  applications  of  phase-field  methodology. 

Key  words:  Phase- field,  Cahn-Hilliard,  Isogeometric  Analysis,  NURBS,  steady 
state  solutions,  isoperimetric  problem 


Preprint  submitted  to  Elsevier  Science 


11  December  2007 


Form  Approved 
OMB  No.  0704-0188 


Report  Documentation  Page 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

DEC  2007 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2007  to  00-00-2007 

4.  TITLE  AND  SUBTITLE 

Isogeometric  Analysis  of  the  Cahn-Hilliard  phase-field  model 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Texas  at  Austin, Institute  for  Computational  Engineering 
and  Sciences, 1  University  Station, Austin, TX, 78712 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

The  Cahn-Hilliard  equation  involves  fourth-order  spatial  derivatives.  Finite  element  solutions  are  not 
common  because  primal  variational  formulations  of  fourth-order  operators  are  well  defined  and  integrable 
only  if  the  finite  element  basis  functions  are  piecewise  smooth  and  globally  Cl-continuous.  There  are  a  very 
limited  number  of  two-dimensional  finite  elements  possessing  Cl-continuity  applicable  to  complex 
geometries,  but  none  in  three-dimensions.  We  propose  Isogeometric  Analysis  as  a  technology  that  possesses 
a  unique  combination  of  attributes  for  complex  problems  involving  higher-order  differential  operators, 
namely,  higher-order  accuracy,  robustness,  two-  and  three-dimensional  geometric  exibility,  compact 
support,  and,  most  importantly,  the  possibility  of  Cl  and  higher-order  continuity.  A  NURBS-based 
variational  formulation  for  the  Cahn-Hilliard  equation  was  tested  on  two-  and  three-dimensional 
problems.  We  present  steady  state  solutions  in  two-dimensions  and,  for  the  first  time,  in  three-dimensions. 
To  achieve  these  results  an  adaptive  time-stepping  method  is  introduced.  We  also  present  a  technique  for 
desensitizing  calculations  to  dependence  on  mesh  refinement.  This  enables  the  calculation  of  topologically 
correct  solutions  on  coarse  meshes,  opening  the  way  to  practical  engineering  applications  of  phase-field 
methodology. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

46 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1  Introduction 


1.1  Phase  transition  phenomena:  the  phase- field  approach 


Two  different  approaches  have  been  used  to  describe  phase  transition  phenom¬ 
ena:  sharp-interface  models  and  phase-field  (diffuse-interface)  models.  Tradi¬ 
tionally,  the  evolution  of  interfaces,  such  as  the  liquid-solid  interface,  has  been 
modeled  using  sharp-interface  models  [37,70].  This  entails  the  resolution  of  a 
moving  boundary  problem.  Thus,  the  partial  differential  equations  that  hold 
in  each  phase  (for  instance,  describing  mass  conservation  and  heat  diffusion) 
have  to  be  solved.  These  equations  are  coupled  by  boundary  conditions  on 
the  interface,  such  as  the  Stefan  condition  demanding  energy  balance  and  the 
Gibbs-Thomson  equation  [6,21,59].  Across  the  sharp  interface,  certain  quan¬ 
tities  (e.g.,  the  heat  flux,  the  concentration  or  the  energy)  may  suffer  jump 
discontinuities.  The  free-boundary  (sharp-interface)  description  has  been  a 
successful  model  in  a  wide  range  of  situations,  but  it  also  presents  complica¬ 
tions  from  the  physical  [2]  and  computational  [15]  points  of  view. 

Phase-field  models  provide  an  alternative  description  for  phase-transition  phe¬ 
nomena.  The  phase- field  method  has  been  used  to  model  foams  [32] ,  describe 
solidification  processes  [10,63],  dendritic  flow  [49,53],  microstructure  evolution 
in  solids  [34],  and  liquid- liquid  interfaces  [61].  For  recent  reviews  of  phase-field 
methods  the  reader  is  referred  to  [16,30]. 

The  key  idea  in  phase-field  models  is  to  replace  sharp  interfaces  by  thin  tran¬ 
sition  regions  where  the  interfacial  forces  are  smoothly  distributed.  Explicit 
front  tracking  is  avoided  by  using  smooth  continuous  variables  locating  grains 
or  phase  boundaries. 

Phase-field  models  can  be  derived  from  classical  irreversible  thermodynam¬ 
ics  [40].  Utilizing  asymptotic  expansions  for  vanishing  interface  thickness,  it 
can  be  shown  that  classical  sharp- interface  models,  including  physical  laws  at 
interfaces  and  multiple  junctions,  are  recovered  [36,38].  In  order  to  capture 
the  physics  of  the  problem,  the  transition  regions  (diffuse  interfaces)  in  the 
phase-field  models  have  to  be  extremely  thin. 

The  use  of  diffuse-interface  models  to  describe  interfacial  phenomena  dates 
back  to  Korteweg  [55]  (1901),  Cahn  and  Hilliard  [12]  (1958),  Landau  and 
Ginzburg  [56]  (1965)  and  van  der  Waals  [72]  (1979). 

The  Cahn-Hilliard  phase-field  model  is  normally  used  to  simulate  phase  seg- 
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regation  of  a  binary  alloy  system,  but  many  other  applications,  such  as,  image 
processing  [23],  planet  formation  [69]  and  cancer  growth  [33]  are  encountered 
in  the  literature. 


1.2  Numerical  methods  for  the  Cahn- Hilliard  phase-field  model 
1.2.1  Spatial  discretization 

The  Cahn-Hilliard  equation  involves  fourth-order  spatial  partial-differential 
operators.  Traditional  numerical  methodologies  for  dealing  with  higher-order 
operators  on  very  simple  geometries  include  finite  differences  (see  applica¬ 
tions  to  the  Cahn-Hilliard  equation  in  [35,67])  and  spectral  approximations 
(solutions  to  the  Cahn-Hilliard  equation  can  be  found  in  [58,60,76,77]).  In 
real-world  engineering  applications,  simple  geometries  are  not  very  relevant, 
and  therefore  more  geometrically  flexible  technologies  need  to  be  utilized.  It 
is  primarily  this  reason  that  has  led  to  the  finite  element  method  being  the 
most  widely  used  methodology  in  engineering  analysis.  The  primary  strength 
of  finite  element  methods  has  been  in  the  realm  of  second-order  spatial  op¬ 
erators.  The  reason  for  this  is  variational  formulations  of  second-order  opera¬ 
tors  involve  integration  of  products  of  first-derivatives.  These  are  well  defined 
and  integrable  if  the  finite  element  basis  functions  are  piecewise  smooth  and 
globally  C°-continuous,  which  is  precisely  the  case  for  standard  finite  element 
functions.  On  the  other  hand,  fourth-order  operators  necessitate  basis  func¬ 
tions  that  are  piecewise  smooth  and  globally  C1 -continuous.  There  are  a  very 
limited  number  of  two-dimensional  finite  elements  possessing  ^-continuity 
applicable  to  complex  geometries,  but  none  in  three-dimensions  (see  [66]  for  a 
recent  study  on  C 1  -continuous  finite  elements).  As  a  result,  a  number  of  differ¬ 
ent  procedures  have  been  employed  over  the  years  to  deal  with  higher-order 
operators.  All  represent  theoretical  and  computational  complexities  of  one  de¬ 
gree  or  another.  Unfortunately,  it  may  be  said  that  after  50  years  of  finite 
element  research,  no  general,  elegant  and  efficient  solution  of  the  higher-order 
operator  problem  exists. 

For  the  above  reasons,  finite  element  solutions  to  the  Cahn-Hilliard  equation 
are  not  common.  The  most  common  way  to  solve  this  equation  in  finite  ele¬ 
ment  analysis  has  been  with  mixed  methods  [4,5,9,24-26,31]  rather  than  the 
use  of  C 1  -continuous  function  spaces  [28].  Recently,  a  discontinuous  Galerkin 
(DG)  formulation  has  been  proposed  (see  [75]).  All  of  these  methods  lead  to 
the  introduction  of  extra  degrees  of  freedom  in  addition  to  primal  unknowns. 
As  a  consequence,  an  alternative  approach  is  desirable.  Perhaps,  the  most  ef¬ 
ficient  procedure  developed  to  date  is  the  so-called  continuous/discontinuous 
Galerkin  (CDG)  method  [29,74],  In  this  method,  standard  C°-continuous  finite 
element  basis  functions  are  used  in  conjunction  with  a  variational  formulation 
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that  maintains  (^-continuity  weakly  through  use  of  discontinuous  Galerkin 
operators  on  derivatives.  This  eliminates  extra  degrees  of  freedom  at  the  price 
of  the  inclusion  of  the  discontinuous  Galerkin  operators  which  change  the  data 
structure  from  the  normal  one  based  solely  on  element  interior  contributions  to 
one  in  which  element  edge  or  surface  contributions  are  additionally  required. 
Nevertheless,  due  to  the  reduction  in  degrees  of  freedom,  this  method  seems 
to  have  the  advantage  over  others  previously  proposed. 


Recently,  a  new  methodology,  Isogeometric  Analysis,  has  been  introduced  that 
is  based  on  recent  developments  in  computational  geometry  and  computer 
aided  design  (CAD)  [46].  Isogeometric  analysis  is  a  generalization  of  finite  el¬ 
ement  analysis  possessing  several  advantages:  1)  It  enables  precise  geometric 
definition  of  complex  engineering  designs  thus  reducing  errors  caused  by  low- 
order,  faceted  geometric  approximation  of  finite  elements.  2)  It  simplifies  mesh 
refinement  because  even  the  coarsest  model  precisely  represents  the  geometry. 
Thus,  no  link  is  necessary  to  the  CAD  geometry  in  order  to  refine  the  mesh, 
in  contrast  with  the  finite  element  method,  in  which  each  mesh  represents  a 
different  approximation  of  the  geometry.  3)  It  holds  promise  to  simplify  the 
mesh  generation  process,  currently  the  most  significant  component  of  analysis 
model  generation,  and  a  major  bottleneck  in  the  overall  engineering  process. 
4)  The  /c-refinement  process,  unique  to  isogeometric  analysis  among  geomet¬ 
rically  flexible  methodologies,  has  been  shown  to  possess  significant  accuracy 
and  robustness  properties,  compared  with  the  usual  /^refinement  procedure 
utilized  in  finite  element  methods  [3,20]. 


^-refinement  is  a  procedure  in  which  the  order  of  approximation  is  increased, 
as  in  the  p- method,  but  continuity  (i.e.,  smoothness)  is  likewise  increased, 
in  contrast  to  the  p-method.  Isogeometric  analysis,  thus,  presents  a  unique 
combination  of  attributes  that  can  be  exploited  on  problems  involving  higher- 
order  differential  operators,  namely,  higher-order  accuracy,  robustness,  two- 
and  three-dimensional  geometric  flexibility,  compact  support,  and,  most  im¬ 
portantly,  C1  and  higher-order  continuity.  In  addition,  higher-order  continuity 
is  achieved  without  introducing  extra  degrees  of  freedom.  These  properties 
open  the  way  to  application  to  phase-field  models.  Herein,  we  report  our  initial 
efforts  to  simulate  higher-order  operators  using  isogeometric  analysis.  Higher- 
order  operators  are  encountered  in  biomedical  applications  and  in  many  areas 
of  engineering,  such  as,  for  example,  liquid-liquid  flows,  liquid-vapor  flows, 
emulsification,  cancer  growth,  rotation-free  thin  shell  theory,  strain-gradient 
elastic  and  inelastic  material  models,  and  dynamic  crack  propagation,  etc.  The 
simplicity  of  isogeometric  analysis  compared  with  many  procedures  that  have 
been  published  in  the  literature  is  noteworthy.  We  believe  it  may  prove  an 
effective  procedure  for  solving  problems  of  these  kinds  on  complex  geometries. 
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1.2.2  Time  discretization 


The  time  integration  of  the  Cahn-Hilliard  equation  is  not  trivial.  The  non¬ 
linear  fourth-order  term  imposes  severe  time-step  size  restrictions  for  explicit 
methods,  thus  mandating  the  use  of  implicit  or  (at  least)  semi-implicit  algo¬ 
rithms.  Under  the  non-realistic  hypothesis  that  assumes  the  mobility  to  be 
constant,  the  fourth-order  term  of  the  Cahn-Hilliard  equation  becomes  linear. 
One  can  take  advantage  of  this  fact  by  using  a  semi-implicit  time  integrator 
that  treats  the  fourth-order  term  implicitly,  while  the  non-linear  second-order 
term  is  treated  explicitly  [77,43] .  This  technique  allows  a  somewhat  larger  time 
step  than  explicit  methods  while  avoiding  the  use  of  nonlinear  solvers.  How¬ 
ever,  in  this  paper  we  are  interested  in  the  thermodynamically  relevant  case, 
where  the  mobility  depends  on  concentration  and  the  fourth-order  term  is  no 
longer  linear.  As  a  consequence,  we  will  use  a  fully  implicit  time  integration 
scheme  which  requires  the  use  of  a  nonlinear  solver. 

Adaptive  time  stepping  is  of  prime  importance  because  the  dynamic  response 
of  the  Cahn-Hilliard  equation  intermittently  experiences  fast  variations  in 
time.  The  usual  approach  presented  in  the  literature  for  simplified  versions 
of  the  Cahn-Hilliard  equation  has  been  to  use  a  few  (2  or  3)  different  time- 
step  sizes  during  the  simulation  [14].  These  time  steps  are  not  selected  by 
means  of  accuracy  criteria,  but  by  using  approximate  theories  of  the  late¬ 
time  behavior  of  the  Cahn-Hilliard  equation  [73] .  In  this  paper  we  propose  an 
adaptive-in-time  method  where  the  time  step  is  selected  by  using  an  accuracy 
criterion.  This  allows  us  to  reduce  the  compute  time  by  factors  of  hundreds 
while  ensuring  that  sufficient  time  accuracy  is  achieved.  (Another  approach 
that  has  been  used  in  the  literature  to  speed  up  the  solution  is  the  use  of 
multigrid  methods  [51,52],  which  is  not  pursued  in  this  work.) 


2  The  strong  form  of  the  Cahn-Hilliard  equation 


Let  Cl  C  be  an  open  set,  where  d  =  2  or  3.  The  boundary  of  0,  assumed 
sufficiently  smooth,  is  denoted  T.  The  unit  outward  normal  to  T  is  denoted 
n.  We  assume  the  boundary  T  is  composed  of  two  complementary  parts,  T  = 
r9  U  rs.  A  binary  mixture  is  contained  in  0  and  c  denotes  the  concentration 
of  one  of  its  components.  The  evolution  of  the  mixture  is  assumed  governed 
by  the  Cahn-Hilliard  equation.  In  strong  form,  the  problem  can  be  stated  as: 
find  c:!1x(0,T)hR  such  that 
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Be 

-  =  V  •  (MCV(/XC  -  AAc)) 

in 

ft  x  (0, T), 

(i.i) 

c  =  g 

on 

r*x(o,T), 

(1.2) 

McV(/ic  —  AAc)  •  n  =  s 

on 

rs  x  (o,  t), 

(1.3) 

McXVc  ■  n  =  0 

on 

r  x  (o,  t), 

(1.4) 

c{x,t )  =  Cq(x) 

in 

ft. 

(1.5) 

where  Mc  is  the  mobility,  //c  represents  the  chemical  potential  of  a  regular 
solution  in  the  absence  of  phase  interfaces  and  A  is  a  positive  constant  such 
that  a/A  represents  a  length  scale  of  the  problem.  This  length  scale  is  related 
to  the  thickness  of  the  interfaces  that  represent  the  transition  between  the  two 
phases. 

Remarks: 

(1)  In  most  of  the  existing  analytic  studies,  as  well  as  numerical  simulations, 
the  mobility  is  assumed  to  be  constant.  However,  according  to  thermo¬ 
dynamics  [11],  it  should  depend  on  the  mixture  composition.  This  depen¬ 
dence  might  produce  quite  important  changes  of  the  coarsening  kinetics. 
In  this  paper  we  consider  the  commonly  adopted  relationship 

Mc  =  Dc{l-  c )  (2) 

in  which  D  is  a  positive  constant  which  has  dimensions  of  diffusivity,  that 
is,  length2/time.  This  relationship  appeared  in  the  original  derivation  of 
the  Cahn-Hilliard  equation  [11]  and  is  commonly  referred  to  as  degenerate 
mobility,  as  pure  phases  have  no  mobility.  We  observe  that  relation  (2) 
restricts  the  diffusion  process  primarily  to  the  interfacial  zones,  which 
is  precisely  what  happens  in  many  physical  situations  where  movement 
of  atoms  is  confined  to  the  interfacial  region  [11],  The  reader  is  referred 
to  the  paper  by  Elliott  and  Garcke  [27]  for  a  proof  of  the  existence  of 
weak  solutions  of  the  Cahn-Hilliard  equation  with  degenerate  mobility. 
Further  information  about  the  regularity  of  the  solutions  can  be  found  in 
[50].  Numerical  simulations  of  the  Cahn-Hilliard  equation  with  degenerate 
mobility  are  reported  on  in  [4,74,75]. 

(2)  The  function  pc  is  a  highly  nonlinear  function  of  the  concentration  rep¬ 
resenting  the  chemical  potential  of  a  uniform  solution  [12].  It  is  usually 
approximated  by  a  polynomial  of  degree  three.  In  this  paper  we  consider 
the  thermodynamically  consistent  function,  namely 

^  =  'gjlog7V7  + 1  ~2c  (3> 

where  6  =  Tc/T  is  a  dimensionless  number  which  represents  the  ratio 
between  the  critical  temperature  Tc  (the  temperature  at  which  the  two 
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phases  attain  the  same  composition)  and  the  absolute  temperature  T. 


2. 1  The  energy  functional  for  the  Cahn- Hilliard  equation 


An  important  feature  of  the  Cahn-Hilliard  model  is  the  existence  of  an  energy 
functional  given  by  the  Ginzburg-Landau  free  energy,  namely 

£{c)=  [  (Tc  +  Ts)da;  (4) 

where  Tc  is  the  chemical  free  energy  and  a  surface  free  energy  term.  Ac¬ 
cording  to  the  original  model  of  Cahn  and  Hilliard  [12,13],  the  surface  free 
energy  is  given  by 

=  AC;-A||Vc||2  (5) 

while  the  chemical  free  energy  takes  the  form 

H/c  =  NkT(c\ogc  +  (1  —  c)  log(l  —  c))  +  Nuc(  1  —  c)  (6) 

where  N  is  the  number  of  molecules  per  unit  volume,  k  is  Boltzmann’s  constant 
and  uj  is  an  interaction  energy,  which,  for  a  system  with  a  miscibility  gap,  is 
positive  and  is  related  to  the  critical  temperature  by 


u  =  2  kTc  (7) 

For  9  =  T'c/T  >  1,  the  chemical  free  energy  is  non-convex,  with  two  wells, 
which  drive  phase  segregation  into  the  two  binodal  points  (values  of  c  that 
minimize  the  chemical  free  energy).  For  6  <  1  it  has  a  single  well  and  admits 
a  single  phase  only. 

The  chemical  potential  /ic  is  given  by  /ic  =  \kc//(Xu;),  where  T'7  is  the  deriva¬ 
tive  of  Tc  with  respect  to  c.  Due  to  the  complexity  of  the  function  Tc,  some 
simpler  approximations  are  normally  employed.  In  particular,  a  polynomial  of 
degree  four  has  been  used  to  approximate  the  chemical  free  energy  in  most  ana¬ 
lytic  studies  and  numerical  simulations  .  The  paper  by  Debussche  and  Detorri 
[22]  (from  the  analytic  point  of  view)  and  the  papers  by  Wells  et  al.  [74], 
Copetti  et  al.  [19]  and  Xia  et  al.  [75]  (from  the  numerical  point  of  view)  deal 
with  the  issue  of  logarithmic  free  energy.  In  the  present  work  we  will  use  the 
logarithmic  function  given  by  (6). 

Remarks: 

(1)  According  to  the  Cahn-Hilliard  model,  the  concentration  is  driven  to  the 
binodal  points  (those  values  of  c  that  minimize  the  chemical  free  energy) 
and  not  to  the  pure  phases. 
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(2)  The  energy  functional  given  by  (4)  constitutes  a  Lyapunov  functional 
since  some  simple  manipulations  lead  to  the  inequality 


d  E 
dt 


AAc  +  /rc)McV (— AAc  +  He) dec  <  0 


where  Li  is  a  real-valued  function  defined  as  Eft)  =  £(c(-,t)). 


(8) 


2.2  Dimensionless  form  of  the  Cahn- Hilliard  equation 


In  the  numerical  examples  presented  in  this  paper,  we  will  use  a  dimensionless 
form  of  the  Cahn-Hilliard  equation.  To  derive  the  dimensionless  equation,  we 
introduce  non-dimensional  space  and  time  coordinates 

x*  =  x/L0,  t*=t/T0  (9) 

where  L0  is  a  representative  length  scale  and  T0  =  Lf/(D\).  In  the  dimen¬ 
sionless  coordinates,  the  Cahn-Hilliard  equation  becomes 

H  =  V*  ■  (A£V*«  -  A *c))  (10) 

where  M*  =  c(l  —  c)  and  /i*  =  [icL\f\. 

We  will  also  make  use  of  the  dimensionless  Ginzburg- Landau  free  energy  given 
by 


E* 


c )  log(l  —  c)  +  26c{l 


c)  + 


dec* 


(11) 


where  E*  =  E(NkTLl)  1  and 


a  = 


T 2 

zE 

3A 


(12) 


is  a  dimensionless  number  related  to  the  inverse  of  the  thickness  of  the  inter¬ 
faces.  The  thickness  of  the  interface  is  inversely  proportional  to  a1/2. 


Following  [74],  we  will  take  the  value  6  =  3/2  for  the  temperatures  ratio  (this 
corresponds  to  a  physically  relevant  case).  Therefore,  the  value  of  a  completely 
characterizes  our  solutions. 


Remark: 

In  what  follows  we  will  use  the  dimensionless  form  of  the  Cahn-Hilliard  equa¬ 
tion.  For  notational  simplicity,  we  will  omit  the  superscript  stars  henceforth. 


3  Numerical  formulation 


3. 1  Continuous  problem  in  the  weak  form 


We  begin  by  considering  a  weak  form  for  the  Cahn-Hilliard  equation.  Let  V 
denote  the  trial  solution  and  weighting  functions  spaces,  which  are  assumed 
to  be  identical.  At  this  point  we  consider  periodic  boundary  conditions  in 
all  directions.  Therefore,  the  variational  formulation  is  stated  as  follows:  find 
c  e  V  such  that  Vie  €  V, 

B(w,c)  =  0  (13) 

where 

B(w,  c )  =  ^w,  ^  +  (Vw,  Mf\7 pc  +  VMcAc)sj  +  (Aw,  McAc)n  (14) 

being  (•,  -)n  the  £2  inner  product  with  respect  to  the  domain  Q.  The  space 
V  =  Tt2  is  a  Sobolev  space  of  square  integrable  functions  with  square  integrable 
first  and  second  derivatives. 

The  repeated  integration  by  parts  of  equation  (14)  under  the  assumptions 
of  periodic  boundary  conditions  and  sufficient  regularity  leads  to  the  Euler- 
Lagrange  form  of  (14): 


/  r\  v 

(«>,  ^  -  V  •  (MCV (/ic  -  Ac))J  =  0  (15) 

which  implies  the  weak  satisfaction  of  equation  (10). 

We  refer  to  this  formulation  as  the  primal  variational  formulation. 


3.2  The  semidiscrete  formulation 


For  the  space  discretization  of  (13)  we  make  use  of  the  Galerkin  method.  We 
approximate  (13)  by  the  following  variational  problem  over  the  finite  element 
spaces:  find  ch  G  Vh  C  V  such  that  \/wh  EVh  C  V 


B(wh,ch)  =  0 


(16) 


where  wh  and  ch  are  defined  as 
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lr 


nb 

Wh  =  Y1  waNa ,  (17) 

A=  1 


nb 

Uh  uANA-  (18) 

A=  1 


The  NaR  are  the  basis  functions,  and  rq  is  the  dimension  of  the  discrete  space. 
Note  that  the  condition  V^CV  mandates  our  discrete  space  to  be  at  least  H2- 
conforming.  This  requirement  is  satisfied  by  a  NURBS  basis  of  (T-continuity  or 
higher.  In  this  paper  we  consider  rectangular  geometries.  In  this  setting,  three- 
dimensional  NURBS  reduce  to  simple  B-splines  in  the  usual  tensor-product 
format  [46].  An  illustration  of  quadratic  B-spline  basis  functions  for  an  eight 
element  mesh  in  one  dimension  is  presented  in  Figure  1.  The  functions  are  C1- 
continuous  at  knots  and  are  C°°-continuous  elsewhere.  Note  that  the  functions 
are  non-interpolatory  at  knots.  As  a  result,  the  solution  coefficients  in  (18), 
referred  to  as  control  variables ,  are  not  associated  with  the  function  value  at 
nodes,  as  in  conventional  finite  element  analysis.  In  the  variational  methods 
literature  they  are  sometimes  referred  to  as  generalized  coordinates. 


3.3  Time  discretization  and  numerical  implementation 


We  integrate  in  time  using  the  generalized- a  method.  This  method  was  origi¬ 
nally  derived  in  [18]  for  the  equations  of  structural  dynamics  and  subsequently 
applied  to  turbulence  computations  in  [3,7,48].  In  addition,  we  propose  a  time- 
step  size  predictor  algorithm  that  allows  us  to  compute  three-dimensional  sta¬ 
tionary  solutions  in  an  affordable  compute  time. 
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3.3.1  Time-stepping  scheme 

Let  C  and  C  denote  the  vector  of  degrees  of  freedom  of  concentration  and 
concentration  time  derivative,  respectively.  We  define  the  residual  vector  as 


R={Ra}  (19.1) 

RA  =  B(Na,  ch)  (19.2) 


The  algorithm  can  be  stated  as:  given  Cn,  Cn  and  A tn  =  tn+ 1  —  tn:  find  Cn+ i, 
C'n+li  Cn+atf  SUCh  that 


Cn+l  =  +  A  tnCn  +  7Atn((7ri+i  —  C7n), 

C* n  T  ^m(^n+l 


'n+ct 


C'n+af  T  ^/(C^n+l  C' r> 


(20.1) 

(20.2) 

(20.3) 

(20.4) 


where  A tn  is  the  current  time- step  size  and  am ,  and  7  are  real- valued 
parameters  that  define  the  method.  Parameters  am,  a /  and  7  are  selected 
based  on  considerations  of  accuracy  and  stability.  Taking  am  =  ay  =  7  = 
1,  the  first  order  backward  Euler  method  is  obtained.  Jansen,  Whiting  and 
Hulbert  proved  in  [48]  that,  for  a  linear  model  problem,  second-order  accuracy 
in  time  is  achieved  if 


1 

7  —  -  +  am  —  a.f 


and  unconditional  stability  is  attained  if 


(21) 


am>af>l/2  (22) 

Parameters  am ,  ay  can  be  parameterized  in  terms  of  ,  the  spectral  radius 
of  the  amplification  matrix  as  At  — >  00,  which  controls  high-frequency  dissi¬ 
pation  [45]: 


1  /  3 


Pc 


2  \  1  +  Pc 


af 


1  +  Pc 


(23) 
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Setting  7  according  to  (21),  a  family  of  second-order  accurate  and  uncondi¬ 
tionally  stable  time  integration  schemes  is  defined  in  terms  of  the  parameter 
Poo  G  [0, 1].  The  nonlinear  system  of  equations  (20)  is  solved  by  using  Newton’s 
method,  which  leads  to  a  two-stage  predictor-multicorrector  algorithm. 

Predictor  stage:  Set 


C  n+ 1,(0) 
Cn+1,(0) 


7  ~ 
7 


(24.1) 

(24.2) 


where  subscript  0  on  the  left-hand-side  quantities  denotes  the  iteration  index 
of  the  non-linear  solver.  This  predictor  was  shown  to  be  efficient  for  turbulence 
applications  [3,7,48]. 

Multicorrector  stage:  Repeat  the  following  steps  for  i  =  1,2,...,  imax 

(1)  Evaluate  iterates  at  the  a- levels 


C n+otm,(i)  "t  ®m(Cn+l,(i-l)  (25.1) 

C n+ctf  ,(i)  C71.  -|“  CXf(Cj 1)  Cn).  (25.2) 

(2)  Use  these  o-level  iterates  to  assemble  the  residual  and  the  tangent  matrix 
of  the  linear  system 


K(i)ACn+lt(i )  =  —R(i)  (26) 

Solve  this  linear  system  using  a  preconditioned  GMRES  algorithm  (see 
Saad  and  Shultz  [64])  to  a  specified  tolerance. 

(3)  Use  ACn+pi)  to  update  the  iterates  as 


n+ 1 ,  (i)  C'  n+l,(i—l)  T  n+l,{i)i  (27.1) 

Cn+i,(i)  =  Cn+i^i_i)  +  7AfnACn+1;p).  (27.2) 

This  completes  one  non-linear  iteration. 

The  tangent  matrix  in  equation  (26)  is  given  by 
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I\ 


dR{C  n-\-arn ,  ac, 


'n+ctn 


dR(cn+am,cn+af)dc. 


'  n+otf 


dC. 


n+an 


dC 


71+1 


dR{C  n-\-arn  ?  C^+qi^) 


ac 


ctf'yAt, 


9Cn+oif  dCn+i 

dR(Cn+am ,  Cn+a  ^ ) 


n+a„ 


dC 


n+OLf 


(28) 


where  the  iteration  index  i  has  been  omitted  to  simplify  the  notation. 
Remarks: 

(1)  The  value  p^  =  0.5  has  been  shown  to  be  an  efficient  choice  for  turbu¬ 
lence  computations  [7].  We  adopted  this  value  for  all  the  computations 
presented  in  this  paper. 

(2)  We  used  the  consistent  tangent  matrix  in  our  computations.  Two  to  four 
nonlinear  iterations  are  typically  required  to  reduce  the  nonlinear  resid¬ 
ual  to  10-4  of  its  initial  value  in  a  time  step.  The  solution  of  system  (26) 
to  a  tolerance  of  10-4  requires  normally  30  to  40  GMRES  iterations  us¬ 
ing  a  diagonal  preconditioner.  The  authors  are  currently  working  on  the 
development  of  more  efficient  preconditioners. 


3.3.2  Time-step  size  adaptivity 

We  borrowed  ideas  from  embedded  Runge-Kutta  methods  [8,41,71]  to  develop 
this  algorithm.  We  took  advantage  of  the  fact  that  the  generalized- a  method 
becomes  the  backward  Euler  method  when  am  =  a  f  =  7  =  1.  The  adaptive 
time  step  strategy  is  presented  in  Algorithm  1.  The  formula  we  use  to  update 
the  time-step  size  is 

(t  A 1/2 

F(e,  At)  =  p  f  j  At  (29) 

Our  default  values  for  the  safety  coefficient  p  and  the  tolerance  tol  are  those 
suggested  in  [57],  that  is,  p  =  0.9  and  tol  =  10-3. 

The  adaptive  time  stepping  technique  allows  us  to  reduce  the  compute  time  by 
factors  of  hundreds  compared  to  the  compute  time  keeping  the  time-step  size 
constant.  Moreover,  it  provides  an  estimate  of  the  time  integration  accuracy. 

Remark: 

When  Algorithm  1  is  used,  the  computed  solution  will  be  rejected  and  recom¬ 
puted  if  the  accuracy  criterion  is  not  fulfilled.  Typically,  fewer  than  10%  of 
the  time  steps  are  rejected  using  the  safety  coefficient  p  =  0.9. 
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Algorithm  1  Time  step  adaptive  process 

Given:  Cn ,  Cn  and  A tn 
1:  Compute  C^_l  using  backward  Euler  and  A tn 
2:  Compute  C“+1  using  second-order  generalized- a  and  A tn 
3:  Calculate  en+1  =  ||C“  -  C“+1||/||C“+1|| 

4:  if  en+i  >  tol  then 

5:  Recalculate  time-step  size  A tn  < —  F(en+i,  Atn) 

6:  goto  1 

7:  else 

8:  Update  time-step  size  Afn+i  =  F(en+ i,  A tn) 

9:  continue 

10:  end  if 

4  Numerical  results 


In  this  section  we  investigate  the  performance  of  our  spatial  and  temporal 
discretization  strategies  for  the  general  Cahn-Hilliard  model.  We  limit  our 
studies  to  simple  geometrical  domains  in  an  effort  to  focus  our  attention  on 
the  physical  and  numerical  aspects  of  the  problem.  The  domain  of  the  test 
cases  is  a  box  Q  =  (0,  l)d,  where  d  =  2  or  3.  At  the  computational  domain 
boundary,  periodic  boundary  conditions  are  imposed  in  all  directions.  The 
spatial  discretization  is  comprised  of  quadratic  spline  functions  that  are  C1- 
continuous  at  knots.  We  employ  meshes  that  are  uniform  in  all  directions. 

The  higher-order  and  higher-continuity  spline  basis  functions  allow  the  use 
of  a  Galerkin  technique  which  yields  a  simple  methodology.  The  efficiency, 
accuracy  and  robustness  of  the  methodology  enabled  us  to  obtain  the  following 
results: 

(1)  Three-dimensional  solutions  for  the  general  Cahn-Hilliard  equation. 

There  are  few  published  numerical  results  of  the  general  case  that  we  an¬ 
alyze  in  this  paper.  To  our  knowledge,  only  [74]  and  [75]  report  numerical 
solutions  to  this  model,  but  the  examples  are  limited  to  2D  domains  and 
early  times. 

(2)  Long-time  behavior  of  the  solution  in  three  dimensions. 

The  behavior  of  the  stationary  solutions  in  multidimensions  is  not  well 
understood  [1,14,17].  In  most  of  the  works  reported  on  in  the  literature 
only  the  fast  evolution  of  the  concentration  that  takes  place  at  the  be¬ 
ginning  of  the  segregation  process  is  computed,  so  the  coarsening  process 
is  completely  neglected.  This  is  due  to  the  fact  that  the  time  integration 
of  the  Cahn-Hilliard  equation  is  very  time  consuming  and  an  efficient 
algorithm  is  necessary  to  be  able  to  compute  stationary  solutions. 

We  remark  that  in  our  experience  obtaining  stationary  solutions  to  the 
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Cahn-Hilliard  equation  is  much  more  challenging  than  computing  only 
initial  fast  dynamics. 

(3)  Statistical  studies  of  solutions  with  random  initial  conditions. 

The  most  commonly  used  initial  condition  for  the  Cahn-Hilliard  equation 
is 

c0(x)  =  c  +  r,  (30) 

where  c  is  a  constant  (referred  to  as  the  volume  fraction)  and  r  is  a 
random  variable  with  uniform  distribution.  This  fact  makes  difficult  the 
comparison  of  solutions  in  terms  of  the  Ginzburg-Landau  free  energy  (the 
quantity  that  better  describes  the  behavior  of  the  solutions  to  the  Cahn- 
Hilliard  equation)  because  it  is  dependent  on  the  initial  condition.  As 
a  consequence,  some  statistics  are  necessary  to  compare  numerical  solu¬ 
tions.  The  statistics  we  consider  in  this  work  are  the  statistical  moments 
up  to  order  10.  The  k-th  order  statistical  moment  is  defined  as, 

Mk=  [  (c  -  c)kdfl  (31) 

Jn 

In  the  numerical  examples,  unless  otherwise  specified,  r  is  a  random  vari¬ 
able  with  uniform  distribution  in  [—0.05,0.05]. 

(4)  A  study  of  several  values  of  the  volume  fraction. 

Although  the  Cahn-Hilliard  phase-field  model  was  proposed  five  decades 
ago,  there  is  still  a  lack  of  understanding  of  basic  points.  In  particular,  a 
study  of  simulations  for  several  values  of  the  volume  fraction  is  lacking, 
even  for  the  quartic  chemical  potential  [68] .  Studies  of  this  kind  are  fun¬ 
damental  to  understanding  the  model  because  the  topology  of  solutions 
is  strongly  dependent  on  the  volume  fraction  c. 


4-1  Numerical  examples  in  two -dimensions 


In  this  section  we  point  out  the  main  difficulties  involved  in  computing  solu¬ 
tions  to  the  Cahn-Hilliard  equation.  We  pay  special  attention  to  the  computa¬ 
tion  of  stationary  solutions,  which  in  our  experience  is  a  much  more  challeng¬ 
ing  problem  than  the  computation  of  transient  solutions  at  early  and  medium 
times.  We  study  the  following  issues: 

(1)  Convergence  of  the  numerical  solution  under  /^-refinement 

(2)  Dependence  of  the  solution  on  the  volume  fraction  c 
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4-1.1  Convergence  of  the  numerical  solution  under  h-refinement 

We  present  two  test  cases  which  are  defined  by  the  sharpness  parameter  a  and 
the  volume  fraction  c.  We  take  the  value  c  =  0.63  in  both  cases.  The  sharpness 
parameter  takes  the  values  3000  and  6000. 

The  initial  condition  is  generated  using  equation  (30).  To  perform  the  h- 
refinement  study,  the  initial  random  distribution  was  generated  on  the  coarsest 
mesh  and  then  reproduced  on  the  finer  meshes.  Thus,  the  initial  condition  is 
exactly  the  same  on  all  meshes  since  the  solution  spaces  are  nested  when  the 
refinement  is  performed  by  knot  insertion  (see  Hughes,  Cottrell  and  Bazilevs 
[46] ) .  The  reason  for  doing  this  is  that  if  the  initial  condition  was  generated  by 
randomly  perturbing  control  variables  on  each  mesh,  the  Ginzburg-Landau  en¬ 
ergy  would  be  significantly  higher  on  the  finer  meshes,  making  the  comparison 
of  the  energy  on  different  meshes  meaningless. 

We  will  compare  the  solutions  on  the  basis  of  statistical  moments  of  order  2, 

3  and  10  and  the  Ginzburg-Landau  free  energy. 

(a)  a  =  3000 

We  computed  solutions  on  322,  642  and  1282  meshes.  The  finer  meshes  are 
obtained  from  the  coarser  meshes  by  knot  insertion  [20,46].  Figures  2,  3  and 

4  show,  respectively,  the  evolution  of  the  statistical  moments  of  order  2,  3 
and  10.  The  evolution  of  the  free  energy  is  depicted  in  Figure  5.  Figures  2- 

5  show  convergence  under  h-refinement.  That  is,  given  an  initial  condition, 
the  evolution  of  the  statistical  moments  and  energy  converges.  The  maximum 
difference  in  Mw  between  the  642  mesh  and  the  1282  mesh  is  less  than  5%, 
and  the  energy  difference  is  barely  discernible. 

Figure  6  shows  snapshots  of  the  solution  computed  on  a  642  mesh  (the  results 
computed  on  a  1282  mesh  were  indistinguishable  from  the  642  case,  while  the 
solution  computed  on  a  322  mesh  was  significantly  less  accurate).  We  observe 
that  the  concentration  is  driven  to  the  binodal  points  in  a  very  fast  process. 
The  separation  occurs  approximately  between  the  times  t  =  2  •  10-6  and 
t  =  4  •  10-6  and  is  driven  by  the  minimization  of  the  chemical  free  energy 
Tc.  This  explains  the  fast  variation  of  the  Ginzburg-Landau  free  energy  that 
takes  place  in  the  time  interval  [2  •  10_6,4  •  10-6]  (see  figure  5).  After  the 
phase  separation,  the  coarsening  process  starts.  The  coarsening  is  driven  by 
the  minimization  of  the  surface  free  energy  The  representative  time  scales 
of  the  coarsening  process  are  much  larger  than  those  of  the  separation  process. 

Remarks: 

(1)  Numerical  solutions  to  this  example  can  be  found  in  [74]  and  [75]  (the 
initial  condition  is  identical  in  the  statistical  sense).  In  those  references, 
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Fig.  2.  Evolution  of  the  second-order  statistical  moment  of  the  concentration  (M2) 
for  a  =  3000  and  c  =  0.63 


Fig.  3.  Evolution  of  the  third-order  statistical  moment  of  the  concentration  (M3) 
for  a  =  3000  and  c  =  0.63 

the  solution  is  computed  on  finer  meshes  (802  linear  elements  in  [75] 
and  10283  quadratic  elements  in  [74]).  Additionally,  in  [75]  the  Cahn- 
Hilliard  equation  is  reformulated  as  a  five-equation  system  that  results 
in  a  much  larger  number  of  degrees  of  freedom.  In  [74]  is  also  reported 
the  solution  to  this  problem  using  mixed  finite  elements  on  a  802  mesh. 
This  approach  also  requires  the  introduction  of  a  number  of  additional 
degrees  of  freedom  as  compared  to  our  method.  Our  solution  computed 
on  a  mesh  of  642  elements  appears  to  be  at  least  of  equivalent  quality  to 
those  reported  on  [74]  and  [75]. 
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Fig.  4.  Evolution  of  the  tenth-order  statistical  moment  of  the  concentration  (M 10) 
for  a  =  3000  and  c  =  0.63 


Fig.  5.  Evolution  of  the  Ginzburg-Landau  free  energy  (E)  for  a  —  3000  and  c  =  0.63 

(2)  The  curve  that  corresponds  to  the  322 3  mesh  in  Figure  5  is  not  defined  for 
certain  times.  This  is  due  to  the  fact  the  mesh  was  not  fine  enough  for 
this  problem  and  the  numerical  solution  at  those  times  was  outside  the 
physical  range  [0, 1],  the  only  interval  where  the  energy  is  defined.  (The 
authors  anticipate  using  the  Variational  Multiscale  Method  [44,47]  to 
address  the  numerical  approximation  of  problems  where  the  characteristic 
length  scale  of  the  equation  is  clearly  unresolved  by  the  computational 
mesh.) 

(3)  It  is  often  argued  that  it  is  not  possible  to  capture  thin  layers  using  high 


18 


(a)  t  =  2.011  •  1CT6  (b)  t  =  4.243  •  1CT6  (c)  t  =  8.327  •  1CT6 


(d)  t  =  1.605  •  10-5  (e)  t  =  3.118  •  10“5  (f)  t  =  6.098  •  10“5 


(g)  t  =  1.324  •  10  4  (h)  t  =  2.656  •  10  4  (i)  Steady  state 


Fig.  6.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  =  3000,  c  =  0.63.  The  mesh  is  comprised  of  642  quadratic  elements. 

continuity  basis  functions.  In  this  example  a  thin  internal  layer  has  been 
captured  in  an  accurate  and  stable  way  using  C 1  -continuous  basis  func¬ 
tions.  In  Figure  8  we  present  cutlines  of  the  steady  solution  for  constant 
values  of  x  corresponding  to  the  vertical  cutlines  represented  on  the  left 
hand  side  of  Figure  7.  In  Figure  8,  n  and  rn  represent  the  knot  number  in 
the  x  direction  and  y  direction,  respectively  (see  also  Figure  7).  We  sam¬ 
ple  the  solution  at  knots  and  plot  it  using  piecewise  linear  interpolation. 
The  solution  is  monotone  and  the  layer  is  captured  within  4  elements. 
We  also  ran  this  example  using  C3-continuous  quartic  basis  functions.  For 
this  case  we  were  forced  to  use  a  different  initial  condition  because  the 
solution  spaces  are  not  nested  when  A:-refinement  is  performed.  In  Figure 
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n  =  i  n  =  21  n  =  32  n  =  1  n  —  10  n  =  40 


Fig.  7.  Steady  state  solutions  of  the  problem  defined  by  a  =  3000  and  c  =  0.63. 
In  both  pictures  (left  and  right),  the  initial  condition  is  the  same  from  a  statistical 
point  of  view,  but  not  from  a  deterministic  point  of  view.  The  mesh  is  comprised  of 
642  quadratic  C 1  -continuous  elements  for  the  solution  on  the  left  hand  side  and  642 
quartic  C3-continuous  elements  for  the  solution  on  the  right  hand  side.  The  vertical 
lines  on  the  left  and  right  hand  side  pictures  represent  the  cutlines  that  have  been 
plotted  in  Figure  8  and  9,  respectively. 

9  we  represent  cutlines  of  the  steady  state  solution  for  constant  values 
of  x  (this  corresponds  to  the  vertical  cutlines  on  the  right  hand  side  of 
Figure  7). 

(b)  a  =  6000 

This  test  case  is  more  challenging  than  the  previous  one,  because  the  param¬ 
eter  a  is  larger.  This  means  that  the  interfaces  are  thinner,  which,  in  turn, 
requires  a  finer  spatial  mesh.  In  addition,  this  case  is  more  interesting  from 
the  physical  point  of  view,  since  phase-field  models  tend  to  sharp-interface 
models  as  the  thickness  of  the  interfaces  tends  to  zero. 

We  computed  solutions  on  642,  1282  and  2562  meshes.  We  plot  snapshots  of 
the  solution  computed  on  the  642  and  1282  meshes.  The  results  on  the  2562 
mesh  were  indistinguishable  from  those  on  the  1282  mesh. 

In  Figures  10  and  11  we  observe  that  the  solutions  at  early  and  medium  times 
on  the  642  and  1282  meshes  are  very  similar.  Only  at  the  steady  state  do  the  so¬ 
lutions  have  significant  differences.  This  example  shows  the  difficulty  involved 
in  computing  steady-state  solutions  to  the  Cahn-Hilliard  equation.  Obtaining 
stationary  solutions  is  much  more  challenging  than  computing  transient  solu¬ 
tions  at  early  and  medium  times  not  only  from  the  point  of  view  of  the  time 
integration,  but  also  from  the  point  of  view  of  the  spatial  discretization. 

The  modification  of  the  parameter  a  has  not  only  changed  the  thickness  of 
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Fig.  8.  Outlines  of  the  steady  state  solution  using  (^-continuous  basis  functions.  The 
mesh  is  comprised  of  642  quadratic  elements.  The  problem  is  defined  by  a  —  3000 
and  c  —  0.63.  The  solution  is  sampled  at  mesh  knots  and  plotted  using  linear 
interpolation.  Symbols  in  the  plot  correspond  to  element  boundaries.  Note  that 
despite  high  order  and  continuity  of  the  basis  functions,  solutions  are  monotone 
and  the  sharp  layers  are  captured  within  4  elements. 


Fig.  9.  Cutlines  of  the  steady  state  solution  using  C3-continuous  basis  functions. 
The  mesh  is  comprised  of  642  quart ic  elements.  The  problem  is  defined  by  a  =  3000 
and  c  =  0.63.  The  solution  is  sampled  at  mesh  knots  and  plotted  using  linear 
interpolation.  Symbols  in  the  plot  correspond  to  element  boundaries.  Note  that 
despite  high  order  and  continuity  of  the  basis  functions,  solutions  are  monotone 
and  the  sharp  layers  are  captured  within  4-5  elements. 
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(a)  t  =  1.878  •  1CT6  (b)  t  =  4.060  •  10“6  (c)  t  =  8.042  •  10“6 


(d)  t  =  1.598  •  10“5  (e)  t  =  3.270  •  10“5  (f)  t  =  6.274  •  10“5 


(g)  t  =  1.257  •  10  4  (h)  t  =  2.665  •  10  4  (i)  Steady  state 


Fig.  10.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  —  6000,  c  =  0.63.  The  mesh  is  comprised  of  642  quadratic  elements. 

the  interfaces,  but  it  has  also  significantly  modified  the  time  scales  of  the 
problem.  For  example,  the  separation  process  is  much  faster  for  a  =  6000 
than  for  a  =  3000. 


Figures  12,  13  and  14  show  the  evolution  of  the  statistical  moments  of  order 
2,  3  and  10,  respectively.  The  evolution  of  the  free  energy  has  been  depicted 
in  Figure  15.  Figures  12-15  demonstrate  convergence  under  //-refinement. 

Remarks: 

(1)  We  investigated  the  possibility  that  the  parameter  tol,  the  tolerance  used 
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(a)  t  =  1.837  •  1CT6  (b)  t  =  3.821  •  1CT6  (c)  t  =  7.801  •  1CT6 


(d)  t  =  1.641  •  10“5  (e)  t  =  3.246  •  10“5  (f)  t  =  6.190  •  10“5 


(g)  t  =  1.237  •  10  4  (h)  t  =  2.557  •  10  4  (i)  Steady  state 

Fig.  11.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  —  6000,  c  =  0.63.  The  mesh  is  comprised  of  1282  quadratic  elements. 


to  estimate  the  time-step  size,  affected  the  results.  We  ran  cases  with 
tol  =  10-3  (our  default  value)  and  tol  =  10-4.  We  found  no  discernible 
differences  in  the  computed  statistics,  while  the  time-step  size  was  signif¬ 
icantly  smaller. 

(2)  We  analyzed  the  performance  of  our  time-step  size  predictor  by  running 
a  simple  case  taking  a  constant  time-step  size  At  =  10-9  (this  was  less 
than  the  minimal  time-step  size  employed  by  our  time  integrator  for  that 
problem).  We  found  no  discernible  differences  in  the  computed  statistics, 
which  suggests  that  our  adaptive  time-stepping  technique  gives  us  very 
accurate  solutions  at  a  small  fraction  of  the  cost  of  the  constant  time-step 
strategy. 
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Fig.  12.  Evolution  of  the  second-order  statistical  moment  of  the  concentration  (M2) 
for  a  =  6000  and  c  =  0.63 


Fig.  13.  Evolution  of  the  third-order  statistical  moment  of  the  concentration  (M3) 
for  a  =  6000  and  c  =  0.63 

(3)  Our  time-stepping  strategy  enabled  us  to  integrate  the  equations  in  time 
up  to  t  ~  10100,  where  the  solutions  were  considered  steady,  at  a  reason¬ 
able  computational  cost. 

(4)  Most  of  the  time-step  size  selection  criteria  that  have  been  previously 
published  are  based  on  the  free-energy  decay.  The  fulfilment  of  inequality 
(8)  at  the  discrete  level  is  the  condition  more  frequently  asked  for  in  Cahn- 
Hilliard  simulations.  In  our  experience,  this  is  a  very  weak  condition  and 
the  numerical  solution  can  be  completely  wrong  even  if  (8)  is  satisfied  at 
the  discrete  level.  Moreover,  we  found  that  the  evolution  of  the  energy  is 
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Fig.  14.  Evolution  of  the  tenth-order  statistical  moment  of  the  concentration  (M 10) 
for  a  =  6000  and  c  =  0.63 


t 

Fig.  15.  Evolution  of  the  Ginzburg-Landau  free  energy  (E)  for  a  =  6000  and  c  =  0.63 

much  more  dependent  on  the  spatial  resolution  rather  than  the  temporal 
resolution. 

(5)  The  time-step  size  predicted  by  our  method  is  fairly  independent  of  the 
spatial  mesh.  We  plot  in  Figure  16  the  evolution  of  the  time-step  size 
for  the  simulation  of  the  problem  defined  by  a  =  3000  and  c  =  0.63  on 
two  different  meshes.  A  similar  behavior  can  be  observed.  Also,  it  is  an 
instructive  exercise  to  compare  the  evolution  of  the  time-step  size  (Figure 
16)  with  the  evolution  of  the  Ginzburg-Landau  free  energy  (Figure  15). 
We  observe  that  the  separation  process  (identified  in  the  energy  plot  by 
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Fig.  16.  Evolution  of  the  time-step  size  on  a  doubly  logarithmic  plot  for  a  =  3000, 
c  =  0.63 

the  fast  variation  in  the  temporal  interval  [2  •  10-6, 4  •  10-6])  is  computed 
using  very  small  time-step  sizes.  After  the  separation  process  is  finished, 
the  coarsening  process  starts.  The  coarsening  takes  place  at  much  larger 
time  scales  what  is  reflected  by  our  time- step  size  predictor.  Finally,  we 
relate  the  fast  variations  of  the  time-step  size  in  the  last  part  of  the 
simulation,  t  E  [2  •  10-4, 4  •  10-2]  to  the  rough  variations  of  the  free  energy 
that  take  place  at  those  times,  which  we  conjecture  is  due  to  coalescence 
of  bubbles. 

(6)  We  use  Gauss  quadrature  in  all  cases.  We  ran  some  of  the  test  cases  using 
3  and  4  integration  points  per  element  in  each  direction.  The  computed 
statistics  were  indistinguishable. 

(7)  Low-order  statistical  moments  are  not  good  discriminators  in  most  of  the 
cases,  while  higher-order  moments  are,  as  M10  illustrates  in  Figures  4  and 
14.  The  Ginzburg-Landau  free  energy  is  the  most  informative  quantity 
for  comparing  solutions  at  steady  state,  since  the  stationary  solution  can 
be  defined  as  the  one  that  minimizes  the  energy  functional.  However,  the 
energy  at  short  and  medium  times  depends  on  the  initial  condition,  which 
makes  use  of  statistical  moments  necessary. 

(8)  We  ran  examples  with  different  distributions  for  the  random  variable  r. 
The  distributions  were  uniform  in  [-L,  L ],  where  L  =  5-10-2  (our  default 
value),  L  =  5  •  10-4  and  L  =  5  •  10-5.  In  the  first  two  cases,  although 
the  statistics  were  different,  we  found  similar  dynamical  processes.  There 
was  no  phase  separation  for  the  case  L  =  5  •  10-5.  In  our  experience  it  is 
very  important  that  the  pseudo-random  number  generator  is  statistically 
reliable.  We  recommend  that  standard  tests  [54]  be  used  to  verify  its 
quality. 

(9)  We  ran  the  previous  examples  on  the  domain  Q  =  (0,  2)2  to  study  the 


26 


Fig.  17.  Evolution  of  the  second-order  (M2)  and  fourth-order  (M4)  statistical  mo¬ 
ments  of  the  concentration  for  a  =  3000  and  c  =  0.50.  The  mesh  is  comprised  of 
1282  quadratic  elements. 

dependence  of  the  solution  on  the  size  of  the  domain.  The  computed 
statistics  were  very  similar. 


4-1-2  Dependence  of  the  solution  on  c 

We  present  the  numerical  solution  to  the  case  defined  by  c  =  0.50  and  a  = 
3000.  The  mesh  is  comprised  of  1281 2  quadratic  elements.  We  plot  the  evolution 
of  the  statistical  moments  of  order  2  and  4  (Figure  17)  and  3  and  5  (Figure 
18).  The  evolution  of  the  Ginzburg-Landau  free  energy  is  depicted  in  Figure 
19.  We  present  snapshots  of  the  solution  in  Figure  21. 

This  solution  presents  two  main  differences  as  compared  to  those  reported  in 
the  previous  section: 

(1)  The  topology  of  the  solution  is  different  than  for  the  case  c  =  0.63.  This  is 
the  typical  topology  corresponding  to  c  =  0.50.  A  deeply  interconnected 
pattern  is  characteristic  of  this  topology.  In  the  cases  presented  in  Section 
4.1.1  the  masses  of  the  two  phases  were  substantially  different  (c  =  0.63 
and  1  —  c  =  0.37),  which  leads  to  the  phenomenon  of  nucleation.  In  this 
case,  one  finds  irregular  droplets  that  evolve  to  circular  shapes  whose 
characteristic  length  increases  with  time  [39].  Also,  we  observe  that  the 
energy  plot  is  smoother  for  c  =  0.50  than  for  c  =  0.63.  This  is  due  to 
the  fact  that  the  coarsening  process  is  more  continuous  for  this  case,  as 
a  consequence  of  the  absence  of  nucleation. 

(2)  In  the  case  c  =  0.50  the  exact  solution  at  the  steady  state  is  a  strip  instead 
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Fig.  18.  Evolution  of  the  third-order  (M3)  and  fifth  order  (M5)  statistical  moments 
of  the  concentration  for  a  —  3000  and  c  =  0.50.  The  mesh  is  comprised  of  1282 
quadratic  elements. 


Fig.  19.  Evolution  of  the  Ginzburg-Landau  free  energy  for  a  =  3000  and  c  =  0.50. 
The  mesh  is  comprised  of  1282  quadratic  elements. 

of  a  circle.  Geometrical  arguments  support  this.  Further  discussion  about 
this  point  will  be  presented  in  the  sequel. 

In  Figure  20  we  present  the  evolution  in  time  of  the  time- step  size  for  this 
case,  which  is  significantly  different  than  that  in  Figure  16.  The  reason  is, 
again,  the  different  topology  of  the  solution.  In  the  case  defined  by  c  =  0.50 
there  is  no  nucleation  and,  as  a  consequence,  the  coarsening  process  is  more 
continuous,  which  is  reflected  by  our  time-step  size  predictor  (compare  Figure 
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Fig.  20.  Evolution  of  the  time-step  size  on  a  doubly  logarithmic  plot  for  a  =  3000, 
c  =  0.50.  The  mesh  is  comprised  of  1282  quadratic  elements. 

16  with  Figure  20).  This  fact  can  also  be  observed  in  the  plot  of  the  Ginzburg- 
Landau  free  energy  (Figure  19)  which  is  significantly  smoother  than  that  of 
Figure  5.  The  only  sign  of  nucleation  is  observed  in  the  (g)  and  (h)  snapshots 
of  Figure  21,  corresponding  to  local  minima  in  the  time-step  size  plot  (Figure 
20)  and  with  small  but  rapid  variations  in  the  Ginzburg-Landau  free  energy 
(see  Figure  19). 


4-1.3  Topology  of  the  steady-state  solution  and  the  isoperimetric  problem 

The  interest  in  stationary  solutions  to  the  Cahn-Hilliard  equation  (besides  the 
multiple  applications  of  this  model  in  material  science)  is  its  relation  with  the 
periodic  isoperimetric  problem  [42],  which  is  one  of  the  major  open  problems  in 
geometry.  It  was  shown  in  [62,65]  that  minimizers  of  the  Ginzburg-Landau  free 
energy  (under  an  appropriate  rescaling)  converge  to  solutions  of  the  isoperi¬ 
metric  problem  when  a  — >  oo  and  9  — >  oo. 

In  a  2D  periodic  square  (2D  flat  torus)  the  solution  of  the  periodic  isoperimet¬ 
ric  problem  is  well  known:  the  minimizers  are  either  a  circle  (when  0  <  c  <  l/V 
or  1  — 1/7 r  <  c  <  1)  or  a  strip  (when  1/7 r  <  c  <  1  —  1/ 7r)  but  in  3D  the  periodic 
isoperimetric  problem  remains  open  even  in  a  periodic  cube 1 . 

The  2D  solutions  to  the  periodic  isoperimetric  problem  (a  circle  and  a  strip) 
correspond  to  the  solutions  that  we  obtained  for  c  =  0.63  (circle)  and  c  =  0.50 


1  It  has  been  conjectured  that  the  possible  solutions  are  either  a  sphere,  a  cylinder 
or  two  parallel  planes. 
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(a)  t  =  1.971  •  1CT6 


(b)  t  =  3.764  •  1CT6 


(c)  t  =  7.601  •  1(T6 


(d)  t  =  1.609  •  10“5 


(e)  t  =  3.134- 10“5 


(f)  t  =  6.530  •  10“5 


(g)  t  =  1.390  •  10  4  (h )  t  =  2.474  •  10  4  (i)  Steady  state 

Fig.  21.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  =  3000,  c  =  0.50.  The  mesh  is  comprised  of  1282  quadratic  elements. 

(strip).  These  solutions  were  obtained  for  a  =  3000  and  a  =  6000  in  the  case 
of  the  circle  and  for  a  =  3000  in  the  case  of  the  strip.  In  all  cases  0  =  3/2. 

According  to  the  results  quoted  for  the  isoperimetric  problem,  the  case  given 
by  c  =  0.63,  although  close  to  the  limit  value  1  —  1/tt,  corresponds  to  a  strip, 
while  the  computed  solution  for  the  Cahn-Hilliard  equation  is  a  circle.  This  is 
due  to  9  not  being  large  enough.  Increasing  9  drives  the  binodal  points  closer  to 
the  pure  phases.  In  this  case,  the  strip  is  an  energetically  better  solution  than 
the  circle.  We  feel  it  is  important  that  our  solutions  correspond  to  solutions 
of  the  periodic  isoperimetric  problem.  This  link  is  relevant  to  discriminate 
between  physically  acceptable  and  unacceptable  solutions. 
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Remark: 


If  the  Cahn-Hilliard  phase-field  model  is  to  be  used  for  the  approximation 
of  solutions  of  the  isoperimetric  problem,  the  authors  recommend  the  use 
of  constant  mobility  and  quartic  chemical  potential,  which  also  converges  to 
solutions  of  the  isoperimetric  problem  as  a  — >  oo,  but  requires  a  significantly 
smaller  compute  time. 


4-1-4  Mesh-independent  Cahn-Hilliard  phase-field  model 

The  Cahn-Hilliard  phase-field  model  converges,  in  a  thermodynamically  con¬ 
sistent  fashion,  to  its  corresponding  sharp-interface  model  as  y/\  (the  charac¬ 
teristic  length  scale  of  the  model)  tends  to  zero.  In  order  for  the  Cahn-Hilliard 
phase- field  model  to  be  realistic  for  engineering  applications,  A  has  to  be  ex¬ 
tremely  small.  On  the  other  hand,  if  the  computational  mesh  is  not  fine  enough 
to  resolve  the  internal  layers  whose  size  is  defined  by  the  length  scale  V%  non¬ 
physical  solutions  are  obtained. 

To  desensitize  this  mesh  dependence,  we  propose  to  relate  the  characteristic 
length  scale  of  the  continuous  phase-field  model  to  the  characteristic  length 
scale  of  the  computational  mesh.  In  the  ideal  case  we  would  obtain  the  best 
approximation  to  the  sharp-interface  model  for  a  given  mesh.  Also,  we  are 
seeking  a  method  that  preserves  the  topology  of  the  solution  at  the  steady 
state  independently  of  the  mesh  size  while  the  thickness  of  the  interface  is 
enlarged  according  to  the  spatial  resolution. 

Numerical  results  using 

A  =  rh2,  (32) 

where  r  is  a  dimensionless  constant,  have  shown  the  potential  of  this  approach. 
The  value  of  r  has  been  determined  by  means  of  numerical  examples.  It  turned 
out  that  the  maximum  value  of  r  that  can  be  used  retaining  mesh  invariance 
depends  on  the  average  concentration  c.  The  closer  c  is  to  0.5  the  larger  r  has 
to  be  and,  consequently,  the  thicker  the  interface.  We  show  examples  using 
r  =  1  for  c  =  0.63  and  r  =  2.5  for  c  =  0.50.  The  results  are  encouraging. 

In  the  first  example,  c  =  0.63.  In  Figure  22  we  show  the  solutions  on  uni¬ 
form  meshes  comprised  of  322,  642  and  1282  quadratic  elements.  The  initial 
condition  is  generated  by  randomly  perturbing  control  variables  on  the  322 
mesh  and  then  it  is  exactly  reproduced  on  the  finer  meshes  as  in  the  previous 
numerical  examples.  On  the  left  hand  side,  we  plot  the  solution  using  for  all 
meshes  A  =  rl28-2  where  r  —  1.  We  find  a  strong  dependence  of  the  solution 
on  the  mesh  size.  On  the  right  hand  side  we  plot  the  solution  adapting  A  to 
the  resolution  of  the  computational  mesh  through  use  of  equation  (32)  with 
r  =  1.  In  this  case,  the  topology  of  the  numerical  solution  is  independent  of 
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the  mesh  size  and  the  interface  is  captured  on  all  meshes  within  4-5  elements. 


In  the  second  example,  c  =  0.5.  In  Figure  23  we  show  the  solution  on  uniform 
meshes  comprised  of  322,  642  and  1282  quadratic  elements.  The  initial  condi¬ 
tion  is  generated  in  the  same  way  as  in  the  previous  example.  The  solutions  on 
the  left  hand  side  have  been  computed  using  A  =  rl28-2  for  all  meshes  with 
r  =  2.5.  A  strong  dependence  of  the  solution  on  mesh  size  is  observed.  The 
solutions  on  the  right  hand  side  of  Figure  23  have  been  computed  adapting  A 
to  the  resolution  of  the  computational  mesh  using  equation  (32)  and  r  =  2.5. 
The  topology  of  the  numerical  solution  is  again  invariant  with  respect  to  the 
mesh  size  and  the  interface  is  captured  on  all  meshes  within  4-5  elements. 

The  previous  examples  show  the  potential  of  the  proposed  approach  to  suc¬ 
cessfully  deal  with  problems  where  the  characteristic  length  scale  of  the  contin¬ 
uous  phase-field  model  is  unresolved  by  the  computational  mesh.  We  believe 
that  with  this  technique  phase-field  modeling,  which  has  been  used  heretofore 
primarily  in  scientific  studies,  may  become  a  practical  engineering  technology. 


4-2  Numerical  examples  in  three- dimensions 


The  complexity  involved  in  the  approximation  of  a  3D  solution  to  the  Cahn- 
Hilliard  equation  is  much  greater  than  for  the  2D  problem.  The  topology  of 
the  solution  is  much  more  complex  and  it  experiences  significant  changes  as 
time  evolves.  There  seems  to  be  almost  nothing  known  about  the  steady  state 
solutions  in  3D. 

There  are  few  references  reporting  numerical  solutions  to  the  Cahn-Hilliard 
phase-field  model  with  degenerate  mobility  and  logarithmic  free  energy.  To  our 
knowledge,  the  numerical  solutions  presented  in  the  literature  are  limited  to 
the  early  part  of  the  simulation  in  2D  domains.  We  present  herein  stationary 
solutions  in  3D  domains. 


4.2.1  a  =  200,  c  =  0.63 

We  computed  the  solution  on  a  1283  mesh.  We  plot  the  evolution  of  the  statis¬ 
tical  moments  of  order  2,  4  (Figure  24)  and  3  (Figure  25).  The  evolution  of  the 
Ginzburg- Landau  free  energy  is  depicted  in  Figure  26.  We  present  snapshots 
of  isosurfaces  of  the  concentration  in  Figure  27.  We  observe  that  the  ran¬ 
domly  perturbed  constant  concentration  evolves  to  a  complex  interconnected 
pattern.  We  did  not  find  any  sign  of  nucleation  for  this  example,  in  contrast 
with  the  2D  counterpart  of  this  problem  (see  Section  4.1.1).  The  steady  state 
solution  is  a  cylinder,  which  is  one  of  the  conjectured  solutions  for  the  periodic 
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(a)  A  =  r  128  2.  The  mesh  is  322 


(c)  A  =  rl28  2.  The  mesh  is  642 


(e)  A  =  r  128  2.  The  mesh  is  1282 


(b)  A  =  r32  2.  The  mesh  is  322 


(d)  A  =  r64  2.  The  mesh  is  642 


(f)  A  =  rl28  2.  The  mesh  is  1282 


Fig.  22.  Steady  state  solutions  to  the  problem  defined  by  c  =  0.63.  We  show  the 
solutions  on  uniform  meshes  comprised  of  322  (a)-(b),  642  (c)-(d)  and  1282  (e)-(f) 
quadratic  elements.  On  the  left  hand  side  we  plot  the  solution  using  A  =  rl28-2 
for  all  meshes  with  r  —  1.  The  dependence  of  the  solution  on  the  mesh  size  h  is 
apparent.  On  the  right  hand  side  we  plot  the  solutions  adapting  A  to  the  resolution 
of  the  computational  mesh.  The  topology  of  the  solution  is  invariant  with  respect 
to  the  mesh  size.  The  only  difference  in  the  solutions  on  the  right  hand  side  is  the 
thickness  of  the  interface. 
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(a)  A  =  r  128  2.  The  mesh  is  322 


(d)  A  =  r64  2.  The  mesh  is  642 


Fig.  23.  Steady  state  solution  to  the  problem  defined  by  c  =  0.50.  We  show  the 
solutions  on  uniform  meshes  comprised  of  322  (a)-(b),  642  (c)-(d)  and  1282  (e)-(f) 
quadratic  elements.  On  the  left  hand  side  we  plot  the  solution  using  A  =  rl28-2 
for  all  meshes  with  r  =  2.5.  The  dependence  of  the  solution  on  the  mesh  size  h  is 
apparent.  On  the  right  hand  side  we  plot  the  solutions  adapting  A  to  the  resolution 
of  the  computational  mesh.  The  topology  of  the  solution  is  invariant  with  respect 
to  the  mesh  size.  The  only  difference  in  the  solutions  on  the  right  hand  side  is  the 
thickness  of  the  interface. 
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Fig.  24.  Evolution  of  the  second-order  (M2)  and  fourth-order  (M4)  statistical  mo¬ 
ments  of  the  concentration  for  a  =  200  and  c  =  0.63.  The  mesh  is  comprised  of 
1283  quadratic  elements. 


Fig.  25.  Evolution  of  the  third-order  statistical  moment  of  the  concentration  (M3) 
for  a  =  200  and  c  =  0.63.  The  mesh  is  comprised  of  1283  quadratic  elements. 

isoperimetric  problem. 


4.2.2  a  =  600,  c  =  0.75 

In  the  2D  case  there  are  small  topological  differences  between  c  =  0.63  and 
c  =  0.75,  but  in  the  3D  problem  there  are  significant  differences.  In  the  latter 
case,  we  found  that  one  of  the  species  nucleated.  As  before,  the  steady  state 
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Fig.  26.  Evolution  of  the  Ginzburg-Landau  free  energy  (E)  for  a  =  200,  c  =  0.63. 
The  mesh  is  comprised  of  1283  quadratic  elements. 

is  a  cylinder. 

We  computed  the  solution  on  a  1283  mesh.  We  plot  the  evolution  of  the  statis¬ 
tical  moments  of  order  2,  4  (Figures  28)  and  3  (Figure  29).  The  evolution  of  the 
Ginzburg-Landau  free  energy  is  depicted  in  Figure  30.  We  present  snapshots 
of  isosurfaces  of  the  concentration  in  Figure  31. 


5  Conclusions 


We  presented  a  numerical  methodology  for  the  approximation  of  the  general 
Cahn-Hilliard  phase-held  model.  Our  method  is  based  on  isogeometric  analy¬ 
sis,  which  allows  us  to  generate  the  C 1  -continuous  functions  required  to  solve 
the  Cahn-Hilliard  equation  in  a  primal  variational  framework. 

We  introduced  an  adaptive  time-stepping  algorithm  which  proved  to  be  very 
effective  for  this  highly  nonlinear  problem.  We  were  able  to  compute  steady 
state  solutions  in  two  and  three  dimensions.  The  spatial  resolution  has  a  sig¬ 
nificant  effect  on  the  numerical  solution  at  the  steady  state.  Non-physical 
solutions  are  obtained  at  the  steady  state  when  the  spatial  mesh  is  not  fine 
enough,  a  shortcoming  not  previously  emphasized  in  the  literature.  However, 
this  effect  is  not  observed  in  the  early  dynamics  of  the  problem. 

In  the  future  we  envision  combining  NURBS  discretizations  that  are  at  least 
C 1  -continuous  at  patch  level  and  employing  the  CDG  approach  of  [29,74]  at  the 
patch  boundaries.  We  feel  that  this  would  provide  a  good  compromise  between 
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1  A  ~~~  - 

(a)  t  =  1.965  •  10“6  (b)  t  =  1.183  •  10“3 

W 

(c)  t  =  1.514  •  1CT3  (d)  t  =  2.368  •  1CT3 

(e)  t  =  5.240  •  10  3  (f)  Steady  state 


Fig.  27.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  =  200,  c  =  0.63.  The  mesh  is  comprised  of  1283  quadratic  elements. 
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Fig.  28.  Evolution  of  the  second-order  (M2)  and  third-order  (M3)  statistical  mo¬ 
ments  of  the  concentration  for  a  —  600  and  c  =  0.75.  The  mesh  is  comprised  of 
1283  quadratic  elements. 


Fig.  29.  Evolution  of  the  third-order  statistical  moment  of  the  concentration  (M3) 
for  a  =  600  and  c  =  0.75.  The  mesh  is  comprised  of  1283  quadratic  elements. 

geometrical  flexibility  and  computational  efficiency.  We  note,  however,  that 
progress  in  computational  geometry  suggests  that  C 1  -continuous  basis,  at  least 
defined  almost  everywhere,  may  be  a  real  possibility  in  the  near  future. 

The  computations  we  performed  demonstrated  that  our  method  is  capable 
of  giving  very  accurate  and  stable  results  even  when  the  solution  possesses 
very  thin  layers  that  evolve  and  propagate  over  the  mesh.  From  the  numeri¬ 
cal  analysis  point  of  view,  it  seems  apparent  that  the  length-scale  parameter 
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Fig.  30.  Evolution  of  the  Ginzburg-Landau  free  energy  (E)  for  a  =  600,  c  =  0.75. 
The  mesh  is  comprised  of  1283  quadratic  elements. 


associated  with  the  thickness  of  interfaces  should  be  redefined  in  terms  of  the 
resolution  of  the  mesh.  Studies  adopting  this  approach  resulted  in  solution 
topology  insensitive  to  mesh  refinement  even  for  very  coarse  discretizations, 
whereas  employing  a  fixed  length-scale  parameter  led  to  unphysical  solutions 
on  coarse  meshes.  Further  studies  of  a  theoretical  nature  investigating  this 
point  should  prove  highly  valuable.  Nevertheless,  this  observation  has  already 
illustrated  the  potential  of  elevating  phase-field  modeling  from  the  realm  of 
purely  scientific  interest  to  a  practical  engineering  level. 


There  is  also  room  for  improvement  in  the  Cahn-Hilliard  model.  We  found 
that  in  the  case  of  nucleation  the  coarsening  process  does  not  take  place  only 
through  bubbles  merging.  Some  of  them  vanish  without  any  contact  with  other 
bubbles.  Conservation  of  the  formulation  creates  a  non-local  “worm  hole” 
effect  in  which  distant  bubbles  enlarge  when  other  non-contiguous  bubbles 
vanish.  We  observed  this  behavior  even  for  fairly  large  values  of  the  thickness 
parameter  ( a  ~  25000). 

In  summary,  this  paper  constitutes  a  first  step  in  the  application  of  isogeomet¬ 
ric  analysis  to  phase-fields  models.  Although  our  initial  efforts  were  focused 
on  the  Cahn-Hilliard  equation,  we  feel  that  the  methodology  presented  herein 
shows  the  way  to  applications  to  other  areas  of  engineering  interest  involv¬ 
ing  higher-order  spatial  operators,  such  as,  for  example,  liquid-vapor  flows, 
rotation-free  thin  shell  theory,  strain-gradient  elastic  and  inelastic  material 
models,  and  dynamic  crack  propagation. 
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(a)  t  =  3.063  •  10“6 


(b)  t  =  1.114  •  10“3 


(c)  t  =  1.236  •  10"3 


(d)  t  =  2.035  •  10“3 


(e)  t  =  4.168  •  10  3  (f)  Steady  state 

Fig.  31.  Evolution  of  the  concentration  from  a  randomly  perturbed  initial  condition 
for  a  =  600,  c  =  0.75.  The  mesh  is  comprised  of  1283  quadratic  elements. 
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